Effective computation of matrix elements between polynomial basis functions

Two methods of evaluating matrix elements of a function in a polynomial basis are considered: the expansion method, where the function is expanded in the basis and the integrals are evaluated analytically, and the numerical method, where the integration is performed directly using numerical quadrature. A reduced grid is proposed for the latter which makes use of the symmetry of the basis. Comparison of the two methods is presented in the context of evaluation of matrix elements in a non-direct product basis. If high accuracy of all matrix elements is required then the expansion method is the best choice. If however the accuracy of high order matrix elements is not important (as in variational ro-vibrational calculations where one is typically interested only in the lowest eigenstates), then the method based on the reduced grid offers sufficient accuracy and is much quicker than the expansion method.  2004 Elsevier B.V. All rights reserved. In many physical applications the problem of effective evaluation of integrals of a multi-dimensional function often arises. For example, the typical situation in the calculation of molecular rotation-vibration spectra is evaluation of matrix elements, which are integrals of the product of two general polynomial functions (basis functions) and an arbitrary function (molecular potential). Depending on the dimensionality and the number of matrix elements such calculations are often computationally very demanding. It is often the case that the properties of the polynomial functions are well known and there are analytical formulae for the integrals of these functions. Ref. [1] considered the case when the basis functions are a direct product of standard polynomials and advocated the expansion method, i.e. expansion of the potential in the basis and subsequent evaluation of the integrals analytically. In this paper, we propose a method of numerical integration on a reduced grid which makes use of the symmetry properties of the basis functions. Considering an example of a non-direct product basis we argue that our method * Corresponding author. E-mail address: [email protected] (M.M. Law). 0010-4655/$ – see front matter  2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2003.12.007 I.N. Kozin et al. / Computer Physics Communications 165 (2005) 10–14 11 of direct numerical integration can be very competitive to the expansion method and may even be the preferable choice. Our starting point is a model for floppy four-atomic molecules as implemented in the new program WAVR4 [2]. The general description of the approach can be found elsewhere [3–6] and more details about the implementation in Refs. [2,7]. Of the six internal coordinates, the three radial coordinates are treated using the discrete variable representation (DVR) [8]. The DVR approximation for the potential reduces six-dimensional integrals to threedimensional. The other three coordinates are the angles (θ1, θ2, φ) which are represented by a non-direct product basis (1) | j l k J K〉 = P̄ |k−K | j (θ1)P̄ k l (θ2) eikφ √ 2π |J K〉, where J and K are the usual rotational quantum numbers, j and l are angular momenta associated with internal coordinates, k is projection of l on the quantization axis, P̄ k l are normalized associated Legendre functions, and |J K〉 are symmetric top eigenfunctions. All the quantum numbers take only non-negative integer values. The most straightforward approach to evaluate integrals is to use numerical integration and the best way to do it is by using a Gaussian quadrature [9,10]. The idea is to replace the integral of a function by the sum of its functional values multiplied by weighting coefficients. The grid points are choosen so that they are the roots of a certain type of polynomial of order N . Then the quadrature is exact for all polynomials up to degree (2N − 1). If we take, for instance, Legendre type polynomials Pl(x) up to order lmax, we need N = lmax + 1 or more quadrature points to compute the integrals